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Abstract 

A numerical scheme is presented for the solution of the compressible Euler equations over unstructured grids in cylin¬ 
drical and spherical coordinates. The proposed scheme is based on a mixed finite volume / finite element approach. 
Numerical simulations are presented for the explosion problem in two spatial dimensions in cylindrical and spheri¬ 
cal coordinates, and the numerical results are compared with the one-dimensional simulations for cylindrically and 
spherically symmetric explosions. 

Keywords: Compressible flows, Cylindrical/Spherical Coordinates, Explosion problem. Finite Element/Volume 
Methods 


1. Introduction 

Some gasdynamics problems exhibit such strong symmetries that could be more convenient to describe these 
problems in a curvilinear reference, i.e. cylindrical or spherical, rather than in a Cartesian one. These are, e.g., astro- 
physical flows. Inertial Confinement Fusion (ICF) applications, sonoluminescence phenomena and nuclear explosions 
[ 1 ]. 

To compute the numerical solution of the compressible flow equations for these kind of flows, an interesting 
possibility is provided by the use of a mixed finite volume (FV) / finite element (FE) approach [2], For example, in 
viscous flows, it is possible to use the FV and the FE to compute the advection and dissipation terms, respectively, 
within the same algorithm. Such a possibility is expected to be of use in the study of the effect of viscosity on e.g. 
the formation of stable shock fronts and on the determination of the onset and dynamics of Richtmyer-Meshkov 
instabilities in cylindrical and spherical implosions [3], 

The combined use of these two different techniques is made possible by the introduction of suitable equivalence 
conditions that relate the FV metrics, i.e. cell volumes and integrated normals, to the FE integrals. Equivalence con¬ 
ditions relating FV and FE schemes have been derived for Cartesian coordinates in two and three spatial dimensions 
[4, 5] and for cylindrical coordinates in axially symmetric two-dimensional problems [6], In both cited references, 
equivalence conditions are obtained by neglecting higher order FE contributions. Subsequently in [7], equivalence 
conditions for the cylindrical coordinates have been derived for the first time without introducing any approximation 
into the FE discrete expression of the divergence operator, and in [8] the differences between the consistent scheme 
and an alternative one violating the the equivalence conditions have been quantified for the case of one-dimensional 
problems in cylindrical and spherical coordinates. In the present paper the consistent formulation is presented in the 
cylindrical and spherical reference systems. Numerical results for the implosion and the explosion problem are also 
provided to demonstrate the correctness of the proposed approach. 

The present paper is structured as follows. In section 2, the scalar equations in cylindrical and spherical coordinates 
are introduced. In sections 3 and 4 the spatial discretizations of the scalar equation by FE and FV are presented in 
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cylindrical and spherical coordinates respectively. Equivalence conditions between FE and FV metrics are also shown. 
In section 6, numerical simulations are presented for the implosion and the explosion problem in the cylindrical 
coordinates on the R-8 and Z-R (axisymmetric) planes and in the spherical coordinates on the r-<p plane. Numerical 
results are also compared to one-dimensional simulations. In section 7 final remarks and comments are given. 

2. Scalar equation in cylindrical and spherical coordinates 

For simplicity we consider first a multidimensional scalar conservation law. The model equation in a three- 
dimensional cylindrical reference reads 


du dfz 
~dt + ~dZ 


+ \-L {RfR) + TM- Q ' 


( 1 ) 


where t is the time, Z, R and 9 are the axial, radial and azimuthal coordinates, respectively, u - u(Z,R,9, t ) is the 
scalar unknown and f 0 (u) - (f/, /r, fg) is the so-called flux function. A more compact form of the above equation is 
obtained by introducing the divergence operator in three-dimensional cylindrical coordinates V 0 - (•) as follows 


d(u) 

dt 


+ V 0 - / 0 (m) = 0. 


( 2 ) 


The model equation in a three-dimensional spherical reference reads 


du 

dt 


J_ d_ 

r 2 dr 


0 rfr ) 


1 d 
r sin 6 d6 


(sin 0/ e ) + 


1 df<p 

rsin# dtp 


= 0 , 


where r. 8 and tf> are the radial, polar and azimuthal coordinates, respectively, u - u(r , 6 , <p, t) is the scalar unknown 
and f°{u) - ( f r ,fe,f<p) is the flux function. In a compact form, using the divergence operator in three-dimensional 
spherical coordinates V - (•), it can be rearranged as 


d(u) 

dt 


+ V -f un = 0. 


(3) 


3. Finite Volume/Element in cylindrical coordinates 

3.1. Node-pair finite element discretization 

The scalar conservation law (2) is now written in a weak form by multiplying it by the radial coordinate R and by 
a suitable Fagrangian test function <pj e V/, c // 1 (Tl). Integrating over the support Q, of ipi gives 

f R^~dn i+ f RipiV® • f 0 (u) dClj = 0, h e N, (4) 

Jn. dt J n . 

where N is the set of all nodes of the triangulation. Note that by multiplying by R. the numerical singularity of 
the cylindrical reference system is formally removed [6], In the following, to simplify the notation, the infinitesimal 
volume dil — RdR d8dZ is not indicated in the integrals. Integrating by parts immediately gives 

f R<Pi^= f Rf 0 -V^,+ f ipif 0 • V 0 R - f Ripif 0 • nf ( 5 ) 

Jn, ot Jo. Jn, Janf 

where dQf = d O, n d O, with dQ, and dil are the boundary of Q, and of the computational domain Q, respectively, 
and where n 0 = nz_Z + n K R + ng9 is the outward normal versor to Q,-. The scalar unknown is now interpolated as 

m(Z, R , 6 , t ) — M/,(Z, R, 9, t ) = ^ Uk{t)ipu(Z, R, 9), 

keN 
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to obtain the Bubnov-Galerkin approximation of the (1), namely 

V ^M 0 = f Rf 0 (u h ) • V 0 (fi + f ^r(u h )-^R- f RVif 0 
tov, dt dn ‘ dan f 

where Nj is the set of shape functions if / ; whose support G/ f overlaps O, of ifj and where 

f Rywk d£l 0 , 

JQ.ik 


(uh)-n 0 , 


(6) 


M 0 = f 

ik 


with O/j = O, n 0/ f . By resorting the so-called flux interpolation technique [9], the flux function /°(m/,) is now 
expanded using the same shape function iff, e V), as follows 


rw=r 


Y u k(m(z, r,&) f k (m(z ’ R - e) - 


\keN 


) keN 


where/?(f) = f 0 (uk(t)), to obtain 


Y m *^f = f R r^n 0 , 

HN, dt ttl, ktkt? Jdn ft 

where 30 fl , = 30, fl 30/ n 30 and Nf is the set of all boundary nodes of O,. Using the following identities 
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demonstrated in Appendix A, the node-pair representation of (6) is find to be 
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where = A/ /\{f} and Nf* = N a \\i }, and the following FE metric quantities have been introduced 
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In the next section, the corresponding FV metrics are derived. 
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3.2. Edge-based finite volume discretization 

The spatially discrete form of the scalar conservation law (2) is now obtained according to the node-centred finite 
volume approach [10]. To this purpose, the integral form of (2) times the radial coordinate R is enforced over a finite 
number of non-overlapping finite volumes C,-, with boundary dC,. Over each control volume C, the cell-averaged 
unknown is introduced as follows 

u(Z,R,G,t) ^ Uj(t) =— I u(Z,R,6,t), 

Vi JCi 

where V, is the volume of the ;-th cell. Splitting the boundary integral into interface and edge contributions 


'■ dr 



( 8 ) 


where is the set of the finite volume Ck sharing a boundary with C„ excluding C, and where dCik - <9C, n dCk + 
0, k + i, is the so-called cell interface. As it is standard practice, the flux vector is assumed to be constant over each 
cell interface. Under this assumption, the domain and boundary contributions read 


fi<r 


f Rn t - 

= f 0 .y 0 

J ik v ik 

with 

v 0 f Rn 0 

and 

JdC ik 

JdC ik 


Jsc ik ' 


f Rf z • 

■ < - ft • 

f R "t 1 
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with 

v 0 d Z f f Rn 0 , 
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JdCf ' 



respectively. If a second-order centred approximation of the fluxes is considered, namely, f 0 = (f 0 + f 0 )/2, the final 
form of the finite volume approximation of ( 2 ) reads. 


V 0 — = 

'■ dr 


Z 

keNi, 


ff+ft 


+rt-v, -rt 
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• i def 

With V; = 


L k - 
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to be compared with the corresponding FE discretization (7). 


3.3. Finite Element/Volume equivalence 

The equivalence conditions relating the above FV metric quantities and the FE ones defined in the previous section 
are now derived. To this purpose, relevant properties of the FE and FV discretizations are briefly recalled. 

Considering FE metric quantities first, from its definition the vector rj 0 , is asymmetric, namely, tj 0 , = - 17 ?. which 

will be referred in the following as property FEM-a. Property FEM-b is obtained by noting that Z/,-c/v, R-SRf = 
0 which gives immediately 

ir=X>:+ff- a®) 

keNi, t 

Property FEM-c stems from the following identity 

3Lf = f Rip,V 0 -x 0 = f RipiX 0 -n 0 - f ifiiX 0 -R- f (11) 

Jn, Jdnf J n, J n, 


where x 0 is the position vector, x 0 = ZZ +RR(0). By substituting the exact expansion x z 
node-pair representation described in the subsection 4.1, Eq. (11) can be written as 


= ZkeN; Xuifik, and applying 


3 L 0 = 


E 

keNj i L 


' **ik 


•x ; + 
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keNt 


■Xf k +x 0 .fi 0 . 


By substituting property FEM-b in the above identity, one finally obtain property FEM-c as 


v.f = z 


keNi, 
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keNt 


■ •X ik ■ 


( 12 ) 
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Considering now FV metric quantities, from the fact that n 0 = —n 0 over <9C«t, property FVM-a reads v 0 . = -v?., 
which corresponds to the conservation property of the scheme. From the Gauss theorem, one also has 



which, from the definition of FV metric quantities, gives property FVM-b as 


r,'= £ 


keN,, 


(13) 


Property FVM-c is obtained by noting that 


3Vf = f RV 0 -x 0 = (f)Rx 0 -n 0 - f 

Jc, JaCi Jc, 


R-x 0 . 


(14) 


The right hand side of (14) is now computed by means of the FV discretization described in subsection 3.2 as 


3Vf = 


£ 

keNn 


■v 0 .-x 0 -L i +x 0 -v 0 

ik i i ii 


which from property FVM-b becomes 


3Vf= 


keNi, 


Vi I JV • 

- '--v 0 

2 lk 


(15) 


Therefore, a FV approximation can be formally obtained from FE metric quantities defined over the same grid 
points by setting (see properties FEM/FVM-a and -b) 


^0 -~*0 


Vi=Li- 

Note that the mass lumping approximation, 

, d u k d Uj 


~ L 0 —- 
,k dr ' dr 


keNj 


has to be introduced in the Eq. (7) for the equivalence conditions to be applicable. By subtracting Eq. (12) to Eq. (15), 
one finally has 


K =a* + E 


■G- £ 


•x ik - 


keNi, 


keNf, 


(16) 


It is remarkable that, differently from the Cartesian case [4, 5], in the cylindrical reference the FV cell is not coin¬ 
cident with the FE lumped mass matrix. Moreover, the shape of the FV cells that guarantees equivalence with FE 
discretization still remains to be determined. 


4. Finite Volume/Element method in spherical coordinates 

4.1. Node-pair finite element discretization 

In the present section the spatial discretization of conservation law in a spherical reference system is presented. 
The scalar conservation law (3) is multiplied by the quantity rsin0 to formally remove the singularity of the 
spherical coordinate system along the zenith direction. The FE node-pair discretization is obtained as done in the 
cylindrical case [11], namely 
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where the following metric quantities have been introduced 


Ml 

def 

1 r sin Otpiipk d£l°. 
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def 
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4.2. Edge-based finite volume discretization 

The FV spatially discrete form of the scalar conservation law (3) 


fr si 

Jciik 
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multiplied by the quantity r sin 6 reads 


where 
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(18) 


4.3. Finite Element/Volume equivalence 

The equivalence conditions relating the FV and FE scheme are obtained like in the cylindrical case. The following 
properties are introduced for the FE metric quantities 
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4 = -4 

(FEM - b) 

V = 2 4 + £ 
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(FEM - c) 
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Jan 3 
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and for the FV metric quantities 
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( 20 ) 


Therefore, a FV approximation can be formally obtained from FE metric quantities defined over the same grid points 
by setting (see properties FEM/FVM-a and -b) 


4 = 4> v°=£°, 

By subtracting (FEM-c) to (FVM-c), one finally has 


V, = L,. 


v i= u i+ Z Z 


•Xl, 


keN, 


keNf. 


( 21 ) 


where, like the cylindrical case, the mass lumping approximation has been introduced. It is clear that also in this case 
he FV cell is not coincident with the FE lumped mass matrix and the shape of the FV cells that guarantees equivalence 
with FE discretization still remains undetermined. 
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5. Fully discrete form of the Euler equations in cylindrical and spherical coordinates 


The Euler equations in cylindrical and spherical coordinates for compressible inviscid flows are now briefly re¬ 
called. The differential form reads 


du 0 

dt 


■ V 0 -f 0 = -s 2 
R 


<9u° 1 

and - +V 0 -f° =- 5 

dt r sin 9 


respectively, where u 0 (Z,R, 9, t) = ( p , m 0 ,E') T and u°(r, 9, <f>, t) = (p , rtf, E') T , with p density, m 0 = (m z , m R , mg ) T , 
itf = (m r , mg, m ( f,) T momentum and E' total energy per unit volume. The flux function and the source term are defined 
as follow in the cylindrical case 


and in the spherical one 
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where fl is the pressure function in terms of the conservative variables. 

The FV spatially-discrete form of the Euler equations in the cylindrical and spherical refrences reads 


du 0 f 0 +f 0 

keNi,± 


and 


( 22 ) 


y e^I = V M 

' dr " *4 2 


ril + TrLi-rr^ + r^ + sfu, 


def Hpf def 

respectively, where the following notation has been used sf - s 0 (\J 0 ), s- — s°( U 0 ), s? - s°(u 0 ). The terms V 0 and 
Vf are computed from the equivalence condition in cylindrical and spherical case, respectively. In the computation, 
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a TVD [12] numerical flux is used, with the van Leer limiter [13]. The fully discrete form of the Euler system is 
obtained by a two-step Backward Differencing Formulae. At each time level, a dual time-stepping technique is used 
to solve the time-implicit problem [14], In all simulation the ideal gas model for the nitrogen (y = c p jc v = 1.39) is 
used. 


6. Numerical results 


In the present section, numerical results for converging and diverging shock waves are presented for both the 
cylindrical and spherical cases. 

The first problem consider is a one dimensional simulation of converging shocks in cylindrical and spherical 
symmetry. This kind of problems can be seen as a special Riemann problem in cylindrical and spherical coordinates 
for which Guderley [15] has showed that the position of the converging shock has a similar solution as follows 


R 

R~c 


A* 



where R — R(t) is the position of the shock as function of the time t, R, is the position of the discontinuity at initial 
time and t c is the time value at which the discontinuity reaches the origin of the domain. The rational number a is the 
exponent of the similar law and A + and A are the coefficients for the converging and diverging shock, respectively. 
All these coefficients can be computed in the case of the polytropic ideal gas [15]. 

The numerical simulation are performed on a one dimensional domain with wall boundary conditions at the two 
boundaries. The initial solution consists of a uniform distribution of the density, the velocity is zero everywhere, while 
the pressure has a discontinuity at the center of the domain with a value of the pressure on the outer part of the domain 
ten times higher that the value on the inner part. 

In figure 1 the Guderley’s law is compared against the numerical results on a grid of 501 (Ar = 2 x 10 3 ) nodes 
with a time step At = 10 4 for both the cylindrical and spherical problems. The numerical and the analytical solutions 
agree very well, even near the origin of the domain where there are strong shocks. After the shocks reflections, there 
are complex interactions between the waves coming to the origin and those reflected that are not taken into account 
by the analytical model. 




(a) (b) 

Figure 1: Position of the converging shock as function of the time for the one-dimensional simulation and for the 
analytical solution of Guderley: (a) cylindrical problem, (b) spherical problem. 





The explosion and the implosion problem are now considered in the case of two spatial dimension. The computa¬ 
tional domain is shown in figure 2, where a representative computational grid is also shown. 



Figure 2: Exemplary grid for the explosion and implosion problems. The grid is the coarse grid. 

Initial conditions for the explosion problem are as follows. The velocity is assumed to be zero everywhere; the 
density is uniform and equal to 1, whereas the pressure is uniform and equal to 10 in a circular region centered at the 
origin with radius r = 0.5. In the remaining portion of the domain, the pressure is uniform and equal to 1. 

The thermodynamics variables are made dimensionless by the corresponding value at rest (outer region), the 
spatial dimensions are made dimensionless by the radius of the domain. Velocity and time are made dimensionless by 
the square root of the reference pressure divided by the reference density; the reference time is the reference length 
divided by the reference velocity. 

In figure 3 and 4 are reported the numerical solution respectively for the cylindrical coordinates system on the 
plane R-8 and spherical coordinates on the plane r-cf >, for the explosion problem. The density isolines at different time 
levels are shown together with the pressure profiles along the axis y = 0. The grid is made with 39 153 nodes 77 587 
triangles and the time step is 2.5 x 10 4 . A shock wave propagates towards the outer boundary of the computational 
domain; the shock wave is followed by a contact discontinuity. A rarefaction wave propagates towards the origin and 
is then reflected outward. Note that the initial corrugation of the shock front, due to the un-even shape of the initial 
discontinuity caused by its discrete representation over an unstructured grid of triangles, is clearly visible also at later 
times. The pressure profiles are compared against reference one-dimensional results for the three different time levels. 
One-dimensional computations were performed over a evenly-spaced grid made of 2 001 nodes, which corresponds 
to an element spacing of 5 x 10 4 . 

A grid dependence study is shown in figure 5. Pressure signals along the y = 0 axis are compared at time t - 0.16 
for three different grid resolutions: the coarse grid is made of 9 551 nodes and 18 745 triangles, the medium one is 
made of 20683 nodes and 40 841 triangles, the fine one is made of 39 153 nodes 77 587 triangles. Numerical results 
are found to be almost independent from the grid resolution. Time step dependence can be appreciated from figure 
6 , where the pressure signals for three different time steps is shown at time t = 0.16 for the medium grid. Numerical 
results are found to be independent from the chosen time step. 

In figure 7 and figure 8, numerical results for the implosion problem are showed in cylindrical and spherical case 
respectively. In the implosion problem the initial conditions are as in the explosion problem where now the high 
pressure region is the outer one and the low pressure region is at r < 0.5. The grid is the fine grid and the time step is 
2.5 x 10~ 4 . A rarefaction wave propagates towards the outer boundary; a shock wave and a contact surface propagates 
inwards. The intensity of the shock increases as it moves towards the origin; when the shock wave is reflected at the 
origin, a region of high pressure/temperature is observed. In the spherical case, this effect is more evident. 
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Figure 3: Density isoline for the cylindrical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline 
corresponds to a density difference of A plp re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: 
the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution. 


Due to the symmetry of the solution and of the computational domain, the spherical implosion and explosion 
problems can be also simulated by axially symmetric Euler equations formulated in a cylindrical coordinate system, 
as done in [7], The spherical and axisymmetric solution correspond to the solution of the same problem on two 
different planes, i.e., the r-(f> plane in the spherical problem and the Z-R plane in the axisymmetric problem, where Z 
is the coordinate along the axis of symmetry and R is the cooridate along the axis normal to the axis of symmetry. The 
solutions in spherical coordinates are shown together with the corresponding solutions in axisymmetric coordinates 
in figure 9 for the implosion problem. In order to reconstruct the complete spherical solution the spherical and 
axisymmetric solutions are represented on each own plane. 
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(c) (d) 

Figure 4: Density isoline for the spherical implosion problem: (a) t = 0.02; (b) t — 0.12; (c) t = 0.16. Each isoline 
corresponds to a density difference of A p/p re f = 0.03. fd) Pressure signal along the y = 0 axis at the same time levels: 
the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution. 


7. Conclusions 

A novel unstructured-grid hybrid finite element/volume method in a cylindrical and spherical reference was pre¬ 
sented. The proposed approach represents an extension of the node-pair scheme earlier for the cartesian coordinates 
and moves from suitable equivalence conditions linking finite element integrals to the corresponding finite volume 
metrics, such as the cell volume or the integrated normals. The equivalence conditions, as done for the first time in 
the case for the cylindrical coordinates [7], were derived here without introducing any approximation and allowed to 
determine all needed finite volume metric quantities from finite element ones. Numerical results are presented for 
one and two spatial dimension compressible flows: these consist in the numerical simulation of the explosion and im- 


11 






Figure 5: Comparison of pressure signals along the y = 0 axis for the explosion problem at time t - 0.16 for different 
grid resolutions against one dimensional simulations. Left: cylindrical case, right: spherical case. 




(a) (b) 

Figure 6: Comparison of pressure signals along the y = 0 axis for the explosion problem at time t — 0.16 for diverse 
time steps against one dimensional simulations. Left: cylindrical case, right: spherical case. 


plosion problems, in which an initial discontinuity in pressure results in the formation of a diverging and converging 
shock, respectively. The computed pressure and density profile agree fairly well with one-dimensional simulation in 
cylindrical and spherical symmetry over a very fine grid. The solutions obtained in spherical reference also agree very 
well with the corresponding solutions in a cylindrical reference where the axisymmetric condition has been used to 
simulate spherical explosions and implosions. 
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Figure 7: Density isoline for the cylindrical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline 
corresponds to a density difference of A plp re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: 
the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution. 


Appendix A. Node-pair finite element for a scalar conservation law 

In this appendix, the following two identities are demonstrated, namely 

Z( f0 , f0 f0 _ f! 

j k + J i _ 0 _ j k ~ J, 

2 ' ” * 2 

.. 

and 
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Figure 8: Density isoline for the spherical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline 
corresponds to a density difference of A p/p re f = 0.03. fd) Pressure signal along the y = 0 axis at the same time levels: 
the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution. 


which together allow to recast the discrete Bubnov-Galekin equation (6) in its node-pair counterpart (7). 

The proof of the identity (A.l) is considered first. The integral on the left hand side of the Eq. (A.l) is assembled 
considering the contributions coming from each element e in the mesh, exploiting the local support property of the 
shape functions, as follows 

X r k ■ f= YjYjf r i- (a.3) 

keNi cefi, keN' 

where £, is the set of the elements having the node i in common and N e is the set of the nodes of element e. The first 
summation on the right-hand side is limited to the elements contained in the support Q, of node i , which are the only 
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(a) 


(b) 



(c) (d) 

Figure 9: Pressure (a), (c) and density (b), (d) contours for the implosion problem at two different time levels: t = 0.1 
first row, t — 0.32 second row. On the horizontal plane is represented the solution of the bidimensional spherical 
problem on the r-cp plane, while on the vertical plane is represented the solution of the cylindrical axisymmetric 
problem on the Z-R plane. 


ones to give a nonzero contribution to integrals containing the function ipi. Note that Q, = \J ee g. ■ Considering now 
the following identity from the Gauss theorem 


which allows to write 


fv 0 (R^<p k )dn 0 = f Ripi<p k n 0 ddtl 0 , 

Jn ' Jan r 

f Ripk^Vt = - f- f ipnp k R + f Ripupkti* 
Jn <• Jn‘ Jn" Jan? 


(A.4) 
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where R = V 'R is the unit vector in the radial direction. Thanks to the previous relation one deduces 

fR<PkV 0 <Pi = l f ^ fR(p k V zl ‘Pi 


^ f Rip k V 0 ipi + ^ (- f RipiV 0 (fik - f <Pi<pkR + f RtPiH 

^ Jo." ^ \ Jn< Jof Jaci' 


-Lrfaz _ _ 
2 lk 2 


<Pi<PkR + 


Rm 

JdQ e 


where in the last equality has been introduced the elemental contributions rf.j^ of the element e to the vector 17 ? 

riT= fR(<PiV*<Pk-<PkV°<Pi), 

such that t} 0 , = 'Z e€ (s i ns k ) ■ By the relation (A.5), the integral (A.3) becomes 

2/*'- fwv.=-i;2/r-|k'+5 fm-i f ««*»')• < 

tti Jo* TZZ.kTZe \ 2 2 2 Jan' I 


On the other hand, from the Eq. (A.5) follows also that 


r - -2 f 

J n 


R(f k V 0 <Pi - wm R + Ripmn 0 , 


which can be recast, using the equation (A.4), as 


i] e : 0 =2 f Rip,V 0 tp k + f <pi<p k R - f Ripiipkti * 
J ci" Jo." Jan' 


Summing up the last relation for all the nodes k belonging to the element Q' and using the fact that YikeW V 0 (fk(x 0 ) = 
0 , Vx 0 e Q e , one obtains 


k(=hJe \ 


<Pi<PkR + Rwmn 0 = 0. 
<? Jan' I 


Summing up the previous relation for all the elements belonging to £, and multiplying by the vector / 0 , follows that 


EE'? + f 

\ Jo? Ja or 


ip k n 0 = 0 . 


Multiplying this relation by 1 /2 and adding it to the right hand side of (A. 6 ) one has 
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By recalling that rf: 0 = 0 for e i (£,- D S k ), that J 7 0 = 2ee(£,n£ t ) and that i] 0 = 0, the right hand side of the last 
equation can be written as 
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hrKt JQ.,I- Vt-Kf 
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that is the relation (A.l). 

Considering now the proof of the identity (A.2), in the left hand side of (A.2) the contribute of the node i is put 
into evidence, namely 


fe/V * feAf , 


Ripupkti 0 + ff 

anf 
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R(pj(p,n 0 . 
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The quantity 
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is now added and subtracted from the right hand side to obtain 

?i<p k n e 
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By recalling that XkcN fk(x 0 ) - l,x 0 e Of, Ve e £, one has 
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which is the relation (A.2). Using the relations (A.l), (A.2) is possible to write (6) as follow 


keN, 
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Zj dr Zj 
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with the metric quantities defined in the subsection 4.1. In the previous relation is necessary to write in node-pair form 
also the second term on the right hand side. By recalling that 

y ft- ffmR = y ft- fpmR + ff- fvmR, 

by adding and subtracting to the previous relation the following quantity 

X y- M 

keN,,* 

one has 

X ft- fpmR = y (ft - ft) • ff’pi-R +y /r • f 

= X ( ff - W ■ fvmR + /f • ff,R (A.9) 

keN,,* Jn ‘ k Jilik 

= Y, ( ff- 0 - \‘PmR + ff-Lf- 
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by substituting this relation into the (A.8) one has the node-pair FE discretization of the scalar equation (2) 


L ^dUi_ __ y 
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keNi,± 
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where the mass matrix has been lumped 


Z „ dll: „ dll: 

Mf — ^ L 0 — , 
d? ' dt 




with Lf - I ] keN . M 0 . 


(A.10) 
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